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I describe a new algorithm for solving nonlinear wave equations. In this ap- 
proach, evolution takes place on characteristic hypersurfaces. The algorithm 
is directly applicable to electromagnetic, Yang-Mills and gravitational fields 
and other systems described by second differential order hyperbolic equations. 
The basic ideas should also be applicable to hydrodynamics. It is an espe- 
cially accurate and efficient way for simulating waves in regions where the 
characteristics are well behaved. A prime application of the algorithm is to 
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C*~> ' Cauchy-characteristic matching, in which this new approach is matched to a 



standard Cauchy evolution to obtain a global solution. In a model problem 
of a nonlinear wave, this proves to be more accurate and efficient than any 
other present method of assigning Cauchy outer boundary conditions. The 
approach was developed to compute the gravitational wave signal produced 
by collisions of two black holes. An application to colliding black holes is 
presented. 
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H ; I. INTRODUCTION 

I will describe a new way to "make waves" . Wave phenomena is treated mathematically 
by hyperbolic equations. Einstein's theory of general relativity is the supreme example. 
It treats gravity as a distortion in the geometry of space-time. Gravitational waves are 
ripples in the geometry which travel through space changing the shape and size of objects 
in their path. At this historic time, Einstein's theory is providing the basis for a new way 
to observe the distant universe. A worldwide network of gravitational wave observatories is 
under construction. It is designed to detect the gravitational waves produced throughout 
the cosmos by collisions between black holes. 

The new way to make waves which I will describe was developed to simulate the pro- 
duction of gravitational waves by using a computer to solve Einstein's equations. But the 
approach applies to any wave phenomena. In order to get the fundamental ideas across, I 
present them in terms of very simple systems. Only near the end will I give in to the urge to 
tell you about the application to black holes. However, even in simple examples, the ideas 
are easiest to understand in terms of space-time language and pictures. 
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The simplest hyperbolic equation is the 1-dimensional advection equation 

(d t + cd x )$ = 0. (1.1) 

It has the general solution $(£, x) = F(x — ct), which describes a wave traveling with velocity 
c with initial waveform F(x) at t — 0. The curves x = xo + ct are called the characteristics 
of the equation. They are the paths in space-time (worldlines) along which you must move 
to "ride the wave" in the sense that $ = const along a characteristic worldline. Since F is 
arbitrary, we can choose F = for x < and F — 1 for x > to obtain a wave with a 
shock front along the characteristic x = ct. In modern parlance, this shock wave carries one 
bit of information. It illustrates the general way information is propagated by waves along 
characteristics in any hyperbolic system. 

The standard way to solve a hyperbolic equation such as Eq. ( [Lip is by means of 
the Cauchy problem. One poses initial Cauchy data $(£ = 0,x) in some spatial region Q 
and then evolves the solution in time. The characteristics causally relate the initial data 
to a unique evolution throughout some domain of dependence D(fl) in space-time. As a 
result, shock fronts, which represent a sudden signal, can only occur across characteristics. 
Hyperbolic equations lead to ordinary differential equations along the characteristics which 
govern the propagation of shock discontinuities. That makes it important for the purpose 
of numerical simulation to enforce the propagation along characteristics as extensively as 
possible. In complicated cases, where the wave velocity has functional dependence c($,x) 
and an analytic solution is not possible, there still exists a well established theory of the 
general properties of hyperbolic systems based upon characteristics |J. 

The method of characteristics is a computational algorithm for evolving Cauchy data 
based upon this theory. As applied to the advection equation (|1.1|), the characteristic veloci- 
ties are calculated at a given discretized time t n = nAt and the algorithm used to determine 
the field at time t n+ i by requiring that $ be constant along the characteristic curve from a 
space-time point P at time level n to point Q at time level n + 1, i.e. 

$ Q = $ P , (1.2) 

where xq = xp + cAt. When the wave velocity is not constant in either space or time, 
the points P and Q along the characteristics do not in general lie exactly on spatial grid 
points, so that interpolations are necessary. But, because there are no sudden changes along 
characteristics (only across them), such interpolations give excellent accuracy by riding the 
wave like a surfer. 

The method of characteristics for the Cauchy problem extends to the generalization of 
Eq. (|1.1|) to any symmetric hyperbolic system of first differential order equations in a multi- 
dimensional space with many evolution variables and with source terms for the creation 
of waves ||. However, there is another classification of hyperbolic systems, which is more 
familiar to physicists, based upon second order differential equations whose principal part 
has the form 

g^dadpQ = 0, (1.3) 

where a is a space-time index, i.e. x a = (t, x, y, z) for 4-dimensional space-time, and a sum 
over the values of the repeated indices is implied. Equation ( |1.3| ) is classified as hyperbolic 
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if the symmetric matrix g a @ has one negative eigenvalue and its remaining eigenvalues are 
positive. Displacements x a — > x a + dx a along characteristic worldlines satisfy 

g a pdx a dx p = 0, (1.4) 

where g a p, called the metric tensor, is the inverse matrix to g a ^. The simplest example is 
the wave equation in one spatial dimension, 




with the displacement along the characteristics satisfying c 2 dt 2 — dx 2 = 0. By rewriting Eq. 
( |1.5|) in the form 

(^dt-d x )(^dt + d s )^ = 0, (1.6) 

a comparison with Eq. ( |1.1| ) shows that the general solution is $(£, x) = F(x — ct) + G(x+ci). 
There are two characteristics x = xo±ct through each spatial point xq. Information is prop- 
agated along these characteristics which now criss-cross the space-time. The conventional 
method of characteristics could be applied to Eq. ( |1.6|) by the standard technique of reduc- 
ing it to a coupled system of first order differential equations. But it is simpler to integrate 
Eq. ( |1.6| ) over a parallelogram £ in space time whose sides are characteristics meeting at the 
corners P, Q, R and S. This gives the 2-dimensional version of Eq. (|1.2|) 

$Q-$ P + $ i? -<l>5 = 0, (1.7) 

as illustrated in Fig. 1 for the case c = 1 where the characteristics x = xq ± t have 45° slope 
in space-time. By using Eq. (|1.7|) , the method of characteristics can be implemented as a 
computational algorithm for the Cauchy evolution of Eq. ( p..5| ) by positioning the point Q at 
time level n + 1, the points P and S at time level n and the point R at time level n — 1. On 
a (t, x) computational grid satisfying Ax = cAt, the characteristics pass through diagonal 
grid points. The Cauchy data, for the wave equation, which consists of the initial values 
<3>(t = 0, x) and d t &(t = 0,x), is used to initialize an iterative evolution scheme at time levels 
n = and n — 1. An exceptional feature is that the one-dimensional wave equation with 
constant wave velocity can be evolved without error by using Eq. ( |1.7| ) as a finite difference 
equation on such a grid. For a variable velocity, the points P, Q, R and S cannot be placed 
exactly on the grid and interpolations are necessary. 

For the wave equation in two or more spatial dimensions, the manner in which character- 
istics determine domains of dependence and lead to propagation equations is qualitatively 
the same. The major difference is that an infinite number of characteristics now pass through 
each point. For the 3-dimensional wave equation, 

[^-d 2 x -d 2 -d 2 }<!> = (1.8) 

the characteristics which pass through the point (x , Vo, ^o) & t time to are straight lines which 
trace out an expanding spherical wavefront of radius c(t — t ) at a given time t; or, from 
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a space-time point of view, these characteristics trace out the 3-dimensional characteristic 
cone 



(x - x o y + (y- y o y + (z- ZoY - c\t - t ) 2 = 0. (1.9) 

In electrodynamics or relativity, the characteristics are light rays and this is called the light 
cone. I will often use that language here. The analogue in hydrodynamics is the Mach cone. 

The future (past) light cone consists of the radially outward (inward) characteristics 
parameterized by t > to (t < to). There is a 2-parameter set of characteristics through 
each point corresponding to the sphere of angular directions (6, <j)) at that point. This leads 
to some arbitrariness in formulating an evolution algorithm for Cauchy data based upon 
the method of characteristics. There are an infinite number of characteristics and associated 
propagation equations which can be used to evolve Cauchy data from time t to time t + At. 

For a practical numerical scheme, it is thus necessary either to average these propagation 
equations appropriately over the sphere of characteristic directions at or select out some finite 
number of characteristics whose propagation equations comprise a nonredundant set. The 
latter approach has been successfully carried out by Butler ||. For the case of plane flow 
of an inviscid fluid (a problem in two spatial dimensions), he formulates an algorithm based 
upon four "preferred" characteristics. In this problem, the geometry is further complicated 
because the characteristics are dynamically dependent upon the fluid variables, in contrast to 



the essentially time independent characteristic cones of Eq. ( |1.9| ). In the numerical scheme, 
the characteristics must themselves be determined by some finite difference approximation. 
That summarizes the Cauchy problem and its solution by the method of characteristics. 



II. CHARACTERISTIC EVOLUTION 

I now present a procedure which avoids the arbitrariness and awkwardness of the method 
of characteristics in a 3-dimensional Cauchy evolution. It is based upon a characteristic 
initial value approach rather than a Cauchy approach. In order to understand the dis- 
tinction, it is essential to view space-time as 4-dimensional, with initial data given on a 
three-dimensional hypersurface. The Cauchy approach is based upon evolution of initial 
data given on a spatial hypersurface in space-time, i.e. points of space at the same instant 
of time t = t , to data at a later time t = t + At. In the characteristic approach, data is 
initially posed on an outgoing characteristic cone (light cone), which defines a hypersurface 
at constant retarded time u = uq. Characteristic evolution then proceeds iteratively to 
characteristic cones u n = uq + nAu. The numerical grid is intrinsically based upon these 
outgoing characteristic hypersurfaces. The algorithm for evolving from retarded time u n to 
u n+ i is based upon characteristics which are uniquely and intrinsically picked out by the 
geometry of the retarded time hypersurfaces. 

This characteristic evolution procedure radically differs from Cauchy evolution. It uses 
concepts developed in the 1960's for studies of general relativity @J^]. These were prompted 
by the inability of the major mathematical tools such as Green's functions and Fourier anal- 
ysis to overcome the difficulties posed by nonlinearity of the equations and the ambiguity of 
the general coordinate freedom in the theory. Its success later motivated a more extensive 
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mathematical treatment of the characteristic initial value problem ||. These new tech- 
niques were especially designed for investigation of gravitational waves. A computational 
algorithm based on this approach has been the most successful in simulating the production 
of gravitational waves from black holes. 
The approach has two novel ingredients: 

• the use of characteristic hypersurfaces to formulate a characteristic initial value prob- 
lem (CIVP) and 

• the use of compactification methods to describe on a finite grid the waves propagating 
to infinity. 

Characteristic hypersurfaces are 3-dimensional sets traced out in space-time by charac- 
teristic worldlines or, equivalently, by their wavefronts. They provide a natural coordinate 
system to describe waves M. It is very fruitful to use an initial value scheme which de- 
scribes time evolution by means of the retarded time coordinate defined by characteristic 
hypersurfaces as a substitute for the familiar Cauchy scheme based upon constant time 
hypersurfaces. As will be shown, this approach uses a completely different form of the 
mathematical equations and the free initial data. 

Compactification methods provide a rigorous description of the radiation field of a source 
observed in the asymptotic limit of going to infinity along a characteristic worldline. ||. 
The key idea is to introduce a new coordinate which ranges over values from to 1 as the 
actual distance from the source ranges from to oo. The hyperbolic equations are rewritten 
in terms of these new coordinates. Asymptotic behavior at the "points at infinity" can 
then be studied in terms of the new coordinate which ranges over finite values. In this 
way, the concept of the radiation zone as an asymptotic limit at infinity is given rigorous 
meaning. Characteristic hypersurfaces are important here since waves travel to infinity 
along characteristic worldlines, not along Cauchy hypersurfaces of constant time. Even for 
field equations as complicated as those of general relativity, this procedure provides a finite 
geometrical description of waves travelling to infinity. The limit points at infinity form a 
boundary to the compactified space-time, which I will refer to as radiative infinity. At a 
given retarded time, this boundary has the topology of a sphere, representing a sphere of 
observers at infinity. 

It should be emphasized that radiative infinity differs drastically from spatial infinity 
(the limit of going to infinity holding time constant). Early considerations of compactifying 
infinite space for computational purposes were discarded because they were based upon 
spatial infinity ||. An attempt to cover infinite space this way by a finite grid at constant 
time fails in a hyperbolic problem because there is necessarily an infinite physical distance 
between a grid point at infinity and its neighbors. This makes it impossible on a spatial 
grid at fixed time to resolve radiation with finite wave length which propagates to infinity. 
In contrast, for a grid constructed on a characteristic hypersurface, the grid points ride the 
wave without noticing its finite wavelength in the approach to radiative infinity. 

In terms of characteristic coordinates u = ct — x and v = ct + x, Eq. (|1.6|) becomes 

d u d v $ = 0, (2.1) 
Here \I/ = d v § satisfies the propagation equation 
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d u V = 



(2.2) 



along the characteristics in the w-direction. This is the essence of how the use of characteristic 
coordinates simplifies the treatment of waves. 

These ideas provide the physical basis for a new computational algorithm [|Kj. I will 



illustrate how it applies to a nonlinear version of the wave equation in 3 spatial dimensions. 
However, this evolution algorithm can be taken over intact to other hyperbolic physical 
systems, including electromagnetic fields, the Yang-Mills gauge fields of elementary particle 
physics, as well as general relativity, because of the common mathematical structure of these 
theories as second differential order hyperbolic equations. Although it has not been explored 
how this approach might be implemented in the case of a first differential order hyperbolic 
system, such as hydrodynamics, the general ideas should be applicable. 



III. CHARACTERISTIC INITIAL VALUE PROBLEM FOR NONLINEAR 

WAVES 

As a simple illustration of the characteristic initial value problem, consider the nonlinear 
scalar wave equation (SWE) in 3-spatial dimensions, which we write in spherical polar 
coordinates (t, r, 9, 0) as 

-dl- d 2 - dl]§ = - ^ r 2 (r$) + ^ = S(9) (3.1) 

r = ^/(x 2 + y 2 + z 2 ), (3.2) 
where c is the wave velocity, L 2 denotes the standard angular momentum operator 

L-» = -*fij£g*>-3». (3.3) 
sin 9 sin 2 9 

and S($>) represents a nonlinear source term. Rather than using ordinary time t, char- 
acteristic evolution uses the "retarded time" coordinate u = ct — r. The outgoing radial 
characteristics are the curves of constant u, 9 and 0. These are the curves in the r-direction, 
holding u = const, shown in the space-time Fig. p]. In (u, r, 9, 0) coordinates, the SWE ( p.l|) 
takes the form 

2d u d r g = d 2 g-^f + rS. (3.4) 

where g = r$. 

The striking feature about Eq.( |3.4|) is that it is only first differential order in retarded 
time u, unlike the more familiar form of the SWE ( |3.1| ) which is of second differential order 
in time t. When data is given on a characteristic initial hypersurface u = Uq, we need only 
specify the initial value of the field and then use Eq. ( [3.4D to calculate its retarded time 
derivative d u Q in order to evolve the initial data. This is in contrast to the conventional 
Cauchy scheme, where both $ and d t <& must be supplied at initial time t = t and Eq. ( |3.1|) 
is the used to compute the second time derivative d 2 Q in order to evolve the data. 
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In a computational implementation of the CIVP, rather than finite differencing Eq. ( |3.4| ) 
directly, it is advantageous to first convert it into an integral equation which is subsequently 
discretized. Using both the outgoing characteristic coordinate u and the ingoing character- 
istic coordinate v = ct + r, the SWE ( |3.4|) takes the form 



T 2 

4d u d v g = --f + rS& (3.5) 

where g = r$. In the u-v plane formed by fixing the angular coordinates (6, 0), we construct 
a parallelogram £ made up of incoming and outgoing radial characteristics which intersect 
at vertices P, Q, R, S as depicted in Fig. [j]. By integrating Eq. ( |3.5| ) over the area £ bounded 
by these vertices, we may establish the identity 
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L 2 g 



9q = 9p + 9s ~ 9r+ 2 / dudr — Z2~ +rS (~) • ( 3 - 6 ) 



This simple identity is the 3-dimensional analogue of Eq. ( |1.7|) , adapted to include a source 
term. It is the starting point for an evolution algorithm which incorporates the essential 
role that characteristics play in the SWE. 

In order to study the far field wave behavior, we transform this equation to the new 
radial coordinate 

x = r/(l + r), < x < 1. (3.7) 

This serves to map an infinite radial domain into a finite coordinate region, and assigns 
infinitely distant radial points to the edge of the coordinate patch (radiative infinity) at 
x = 1, where the radiation signal can be identified. 

A. Numerical Algorithm 

To develop a discrete evolution algorithm, we work on the lattice of points 

u n = nAu (3.8) 
Xt = iAx 
Q,k = (j + ik)A<p 

where a complex stereographic coordinate ( is used to cover the the sphere in two patches 
(North and South) in order to avoid the polar singularities of the spherical coordinates (8, </>). 
We denote the field at these sites by 

9i jk = g(u n ,Xi,(j,k)- (3.9) 

(We will generally suppress the angular indices j and k.) 

With respect to the (u, x) coordinate grid, it is not possible to place the corners P, Q, R 
and S at grid points since the slope of the characteristics in the compactified x coordinate 
depends upon location. As a consequence, the field g at these points must be interpolated 
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from neighboring grid points. The essential feature of the placement of the parallelogram 
on the grid is that the sides formed by the ingoing characteristics intersect adjacent u- 
hypersurfaces at equal but opposite x-displacements from the neighboring grid points, as 
illustrated in Fig. |2|. The field values at the vertices of the parallelogram are obtained by 
quadratic interpolation. Cancellations between the interpolation errors at the four vertices 
yields the accuracy 

9q ~ 9p ~ 9s + 9r = G Q - G P - G s + G R + 0((Ax) 3 Au), (3.10) 

where G represents the exact analytic solution. 

The integral in Eq. fl3.6| ) can be evaluated by treating the integrand as a constant over 
the parallelogram, with value at the center. The radial coordinate of the point at the center 
is r c = (rp + rg)/2. To compute the nonlinear term, the value of g at r c is taken as the 
average g c = (gp + gs)/2, with g P and g$ evaluated from second-order linear interpolations 
over adjacent points on the grid. The angular derivatives in Eq. (|3.6|) are replaced with 
standard second-order-accurate finite difference approximations. L 2 g is calculated on the 
grid points, and the same interpolation procedure is used to obtain the value of L 2 g c . The 
integral term is then approximated by 

L 2 g + r 3 S(g/r)}dudr/r 2 = [-L 2 g c + r 3 S(g c /r c )} / dudr/r 2 

= 2log( r -^)[-L 2 g c + r 3 c S( 9 ^)], (3.11) 
r P r s r c 

where the integrand is accurate to second order in Ax and A<f>. The resulting finite difference 
equation 

9i +1 ( x Q + x P - Xi-2 ~ Xi-i) = 

29i-i( x Q + x P - Xi-2 - Xi) - g?+£(xQ + x P - Xi-i - Xi) 
+{9i+i(x s +xr~ Xi-t - Xi) - 2g?(x 3 + x R - x^ - x i+1 ) 

(x s - Xr) 



+g?_ 1 (x s + x R -Xi- x i+1 )}- 



+ 



drdu 



s 



r 



2 



-2/ \ i JAq/9c 



[Xq - X P ) 

(Ax 



,2 



-L\g c ) + rtSC-) ^ ' . (3.12) 



(Xq - X P 



relates values of g™ +1 with values at neighboring grid points which are either earlier in 
retarded time (g^, <?!+i), or else contemporary but located at smaller radius (g™^-* 9?-i)- 
Consequently, it is possible to move through the grid, computing g™ +1 explicitly by an orderly 
march. This is achieved by starting at the origin at time u n+ \. Field values of g = r$ vanish 
there. Step outward to the next radial point, using Eq. (|3.12| ) for all angular sites on the 
grid, and iterate this march out to radiative infinity thus updating the characteristic cone at 
u n+ i and completing one retarded time step. This march is then iterated in retarded time. 

The algorithm steps g radially outward one cell with a local error of fourth order in grid 
size. This leads to second order global accuracy which is confirmed by convergence tests 



using known analytic solutions [10]. A complete specification of the algorithm would require 



a description of how the startup procedure at the origin is handled and how stereographic 
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coordinates are used to compute angular derivatives in a smooth way. Details of the use of 
North and South stereographic grids are given in Ref. Physical behavior of nonlinear 
waves is treated in Ref. [K|. Construction of an exact nonspherical solution for an S = $ 3 
self-interaction allows calibration of the algorithm in the nonlinear case where physical sin- 
gularities form. The predicted second order accuracy is confirmed right up to the formation 
of the singularity. Other choices of nonlinear potential allow simulation of solitons. 

The Courant-Friedrichs-Lewy (CFL) condition that the numerical domain of dependence 
contain the physical domain of dependence (determined by the characteristics) is a necessary 
condition for convergence of a finite difference algorithm. For a grid point at (u,r,9), there 
are three critical grid points, at (u — Au, r + Ar, 6) and (u — Au, r — Ar, 9 ± AO), which must 
lie inside its past characteristic cone. These gives rise to the inequalities Am < 2Ar and 
Am < — Ar + (Ar 2 + r 2 A0 2 ) 1 / 2 . At large r, the second inequality becomes Am < rA9 and the 
Courant limit on the time step is essentially the same as for a Cauchy evolution algorithm. 
However, near the vertex of the cone, the second inequality gives a stricter condition 

An < KArAO 2 , (3.13) 

where the value of K depends upon the start up procedure at the vertex. For the scalar 



wave equation, these stability limits were confirmed by numerical experiments |T(| and it 
was found that K « 4. 

Local von Neumann stability analysis leads to no constraints on the algorithm. This 
may seem surprising because no analogue of a CFL condition on the time step arises. It 
can be understood in the following vein. The local structure of the code is implicit, since it 
involves 3 points at the upper time level. Implicit algorithms do not necessarily lead to a 
CFL condition. However, the algorithm is explicit in the way that the evolution starts up 
as an outward radial march from the origin. It is this startup procedure that introduces a 
CFL condition. 

Operating within the CFL limit, the algorithm gives a stable, globally second order 
accurate evolution on a compactified grid ||10|| . (In some nonlinear applications, artificial 
dissipation is necessary for stability ||12|| .) Radiative infinity behaves as a perfectly trans- 
mitting boundary so that no radiation is reflected back into the system. Numerical evolution 
satisfies a conservation law relating the loss of energy to the radiation flux at infinity. 

That summarizes the characteristic initial value problem and its implementation as a 
new computational approach to simulate waves. 



IV. CAUCHY-CHARACTERISTIC MATCHING 

Characteristic evolution has many advantages over Cauchy evolution. Its one disadvan- 
tage is caused by the existence of either (i) caustics where neighboring characteristics focus 
or (ii), a milder version of this, cross-over points where two distinct characteristics collide. 
The vertex of the characteristic cone is a highly symmetric example of a point caustic where 
a complete sphere of characteristics focus. I have already discussed how a point focus gives 
rise to a strong limitation imposed by the CFL condition. 

Cauchy-characteristic matching (CCM) is a way to avoid such limitations by combining 
the strong points of characteristic and Cauchy evolution in formulating a global evolu- 
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tion. |7 yi3|jr^| Here I illustrate the application of CCM to the nonlinear wave equation fl3.1|) . 
This problem requires boundary conditions at infinity which ensure that the total energy 
and the energy loss due to radiation are both finite. In a 3-dimensional problem, these 
are the conditions responsible for the proper 1/r asymptotic decay of the radiation fields. 
However, for practical purposes, in the computational treatment of such a system by the 
Cauchy problem, an outer boundary is artificially established at some large but finite dis- 
tance. Some condition is then imposed upon this boundary in an attempt to approximate 
the proper asymptotic behavior at infinity. Such an artificial boundary condition (ABC) 
typically causes partial reflection of the outgoing wave back into the system P]r5|-|r7[ , which 
contaminates the accuracy of the evolution and the radiated signal. Furthermore, nonlinear 
wave equations often display backscattering so that it may not be correct to try to entirely 
eliminate incoming radiation from the numerical solution. The errors introduced by ABC's 
are of an analytic origin, essentially independent of the computational discretization. In 
general, a systematic reduction of the error can only be achieved by simultaneously refining 
the grid and moving the computational boundary to a larger radius, which is computation- 
ally very expensive for three-dimensional simulations. CCM provides a global solution which 
does not introduce error at the analytic level. 

For linear wave problems, a variety of ABC's have been proposed. For recent reviews, see 
Ref's. [fF?|-|20fl. During the last two decades, local ABC's in differential form have been ex- 
tensively employed by several authors [15, |l6l,pl -|25|l with varying success. Some local ABC's 
have been derived for the linear wave equation by considering the asymptotic behavior of 
outgoing solutions |22] ; this approach may be regarded as a generalization of the Sommerfeld 
outgoing radiation condition. Although such ABC's are relatively simple to implement and 
have a low computational cost, their final accuracy is often limited because their simplifying 
assumptions are rarely met in practice [rS,B|. Systematic improvement of the accuracy of 
local ABC's can only be achieved by moving the computational boundary to a larger radius. 

The disadvantages of local ABC's have led to implementation of nonlocal ABC's based 
on integral representations of the infinite domain problem |T8| , p!9|j29|| . Even for problems 
where the Green's function is known and easily computed, such approaches were initially 
dismissed as impractical pi |; however, the rapid increase in computer power has made it 



possible to implement nonlocal ABC's for the linear wave equation even in 3 space dimen- 
sions p6| . For a linear problem, this can yield numerical solutions which converge to the 



exact infinite domain problem as the grid is refined, keeping the artificial boundary at a 
fixed distance. However, due to nonlocality, the computational cost per time step usually 
grows at a higher power of grid size (0(N 4 ) per time step in a 3-dimensional problem with 
0(N 3 ) spatial grid points) than in a local approach [18,19,2^], which is demanding even for 
today's supercomputers. Further, the applicability of current nonlocal ABC's is restricted 
to problems where nonlinearity may be neglected near the grid boundary ||19|| . 

To my knowledge, only a few works have been devoted to the development of ABC's 



for strongly nonlinear problems [18,27,28]. In practice, nonlinear problems are often treated 



by linearizing the governing equations in the far field , using either local or nonlocal linear 



ABC's |19|j20|] . Besides introducing an approximation at the analytical level, this procedure 
requires that the artificial boundary be placed sufficiently far from the strong-field region, 
which sharply increases the computational cost in multidimensional simulations. There 
seems to be no currently available ABC which is able to produce numerical solutions which 
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converge (as the discretization is refined) to the infinite domain exact solution of a strongly 
nonlinear 3- dimensional wave problem, keeping the artificial boundary at a fixed location. 

For such nonlinear problems, CCM produces an accurate solution out to radiative infin- 
ity with effort 0(N 3 ) per time-step. CCM increases the total computational cost only by a 
factor ~ 2 with respect to a pure Cauchy algorithm with a local ABC. The use of numerical 
methods based upon matching a characteristic initial-value formulation and a Cauchy for- 
mulation can effectively remove the above difficulties associated with a finite computational 
boundary. There is no need to truncate space-time at a finite distance from the sources, since 
compactification of the radial coordinate makes it possible to cover the whole space-time 
with a finite grid. In this way, the true radiation zone signal may be computed. Although 
the characteristic formulation has stability limitations in interior region where the character- 
istic hypersurfaces can develop caustics, it proves to be both accurate and computationally 
efficient in the treatment of the exterior, caustic-free region. 

CCM is a new approach to global numerical evolution which is free of error at the analytic 
level. The characteristic algorithm provides the outer boundary condition for the interior 
Cauchy evolution, while the Cauchy algorithm supplies the inner boundary condition for 
the characteristic evolution. Since CCM consists of discretizing an exact analytic treatment 
of the radiation from source to radiative infinity, it generates numerical solutions which 
converge to the exact analytic solution of the radiating system even in the presence of 
strong nonlinearity. Thus, any desired accuracy can be achieved by refining the grid, without 
moving the matching boundary. In practice, the method performs extremely well even at 
moderate resolutions. 

A. CCM for nonlinear waves 

As an illustration of the computational implementation of CCM , consider again the 
nonlinear 3-dimensional wave equation For simplicity, set the velocity c = 1. In 

the standard computational implementation of the Cauchy problem for ( |3.1|) , initial data 
$(t ,x, y,z) and d t $(t , x, y, z) are assigned and evolved in a bounded spatial region, with 
some ABC imposed at the computational boundary. In a characteristic initial-value for- 
mulation the SWE Q3.1|) is reexpressed in the form of Eq. fl3.4|) in terms of g = r$, using 
standard spherical coordinates and a retarded time coordinate u = t — r. The initial data 
g(u ,r,9,4>), on an initial outgoing characteristic cone u = u is then evolved globally out 
to radiative infinity. 

In CCM, is solved in an interior region r < R m using a Cauchy algorithm, while 
a characteristic algorithm solves the retarded coordinate version ( |3.4| ) for r > R m . The 
matching procedures ensure that, in the continuum limit, $ and its gradient are continuous 
across the interface r = R m . This is a requirement for any consistent matching algorithm, 
since a discontinuity in the field or its gradient could act as a spurious boundary source, 
contaminating both the interior and exterior evolutions. 

For technical simplicity, I illustrate the details of the method here for spherically symmet- 
ric waves but it has been successfully implemented in fully nonlinear 3-dimensional problems 
without symmetry. In a 3-dimensional problem without symmetry, the characteristic evolu- 
tion is carried out on an exterior spherical grid, while the Cauchy evolution uses a Cartesian 
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grid covering the interior spherical region. Although a spherical grid could also be used in 
the interior, a Cartesian grid avoids the necessity of cumbersome numerical procedures to 
handle the singularity of spherical coordinates at the origin. A Cartesian discretization in 
the interior and a spherical discretization in the exterior are the coordinates natural to the 
geometries of the two regions. However, this makes the treatment of the interface somewhat 
involved; in particular, guaranteeing the stability of the matching algorithm requires careful 
attention to the details of the inter-grid matching. Nevertheless, there is a reasonably broad 
range of discretization parameters for which CCM is stable |30| . 

With the substitution G = r<fi and the use of spherical coordinates, the spherically 
symmetric version of the SWE ( |3.1| ) reduces to the 1-dimensional wave equation 

d tt G = d rr G + rS. (4.1) 

The initial Cauchy data is G(t ,r) and d t G(t ,r) in the region < r < R m . Together with 
the regularity condition G(t, 0) = 0, these data determine a unique solution in the domain 
of dependence D\- indicated in Fig. The outer boundary of the domain of dependence is 
the ingoing radial characteristic C\- described by r = R m — t + to- The solution cannot be 
constructed throughout the complete interior region r < R m without additional information, 
which can be furnished by giving the value of G on the outgoing characteristic C + described 
by r = R m +t— t Q (see Fig. |3|). In terms of the coordinates u = t—r and r, C + is described by 
u = to — R m In these coordinates, expressing g(u,r) = G(u + r,r), the spherically symmetric 
version of the wave equation (|3.4j ) is 

2d ur g = d„g + rS. (4.2) 



A unique solution of (|4.2j ) is determined by characteristic initial data consisting of the value 
of g on the initial outgoing characteristic Cq+ and on the ingoing characteristic C\-. These 
data determine the solution uniquely throughout the future of C\- and C +, i.e. the region 
D 1+ in Fig. §. 

The matching scheme proceeds as shown in Fig. |3|. First, initial Cauchy data are evolved 
from t to ti throughout the region D t -, which is in its domain of dependence. Next, this 
induces characteristic data on C\- which combined with the initial characteristic data on 
Cq+ allows a characteristic evolution throughout the region bounded in the future by 
the characteristic C\+ . The solution determined from this initial stage induces Cauchy data 
at time t\ in the region r < R m , inside the matching boundary. This process can then be 
iterated to carry out the entire future evolution of the system. 

CCM is a discretized version of this scheme in which the criss-cross pattern of charac- 
teristics inside the radius R m is at the scale of a grid spacing. The discretized evolution 
algorithm consists of the following steps (see Fig. ^) : 

Step 1. Cauchy evolution. The interior integration scheme is implemented on a uniform 
spatial grid = iAr (0 < i < M) with outer radius Rb = MAr. We discretize ( |4.1| ) using 
the standard second-order finite difference scheme 

Gn+l n/^n i /~in—l rin O/^n _j_ (~in 

(Aly " {A^f + r ^ ' [A - 6) 

where G™ = G(t n ,rj), = S(t n ,ri), and t n = t + nAt. The interior evolution is initialized 
by evaluating G° and G\ (0 < i < M) to second order accuracy from the Cauchy initial 
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data. In the nth time step, (fO|) is used to compute 1 < i < M in terms of field 

values G™ -1 and G™. The regularity of $ at r = implies that Gq = for all n. The 
boundary values G\ j[+1 which are required by ( f4.3|) are supplied by the matching procedure 
(step 3). 

Step 2. Characteristic evolution. The characteristic algorithm is implemented on a 
uniform grid based on the dimensionless compactified radial coordinate 

" = 1 -TT7^^ £ " £1 (4 ' 4) 

so that points at radiative infinity (corresponding to n = 1) are included in the grid. In 
order to include one of the technical problems in matching an interior Cartesian grid to an 
exterior spherical grid, let there be a small gap between the outer radius Rb of the Cauchy 
grid and the matching radius R m (which is also the inner radius of the characteristic grid). 
In 3-dimensional Cartesian-spherical matching the outermost Cauchy grid points and the 
innermost characteristic grid points are necessarily distinct. This is represented here by 
a gap R m — Rb = nAr, where k > is an arbitrary parameter. The characteristic grid 
consists of the uniformly spaced points n a = | + a An (0 < a < N v ), where A?? = {2N v )~ l . 
The retarded time levels u = u n for the characteristic evolution are chosen so that u = u n 
intersects the time level t — t n of the Cauchy evolution at the matching radius; therefore, 
u n = t n — R m and An = At. We denote by g% the value of g at n = n Q , n = u n . The initial 
characteristic data consist of < a < N v . 

In the nth iteration of the evolution, we compute the field values at the grid points 
with u = u n using the values of g^ -1 , which are known either from initialization or from 
the previous iteration. As already described, this is done by the characteristic marching 
algorithm based on the integral identity (p.6|). 

At the inner boundary of the characteristic grid (a = 1), the previous scheme must be 
slightly modified, since <?™_ 2 is not defined. For this initial step, PQRS is chosen so that 
Vp = Voi Vq = n u an d gn, gs are approximated by quadratic interpolation in terms of 
<7o _1 , g™ ~ , (?2 _1 5 which have already been computed. Besides these field values, the final 
evaluation of still requires the value of g^ , which is supplied by the matching procedure 
(step 3). 

Step 3. Matching. Numerous schemes are possible in the case of spherically symmetric 
matching. Here I describe one which is stable for a wide range of gap sizes, < k < 2 
and works in the more complicated situation with an interior Cartesian grid and an exterior 
spherical grid. 

The required boundary values G\ I+l and g$ are computed by radial interpolations at 
constant t, using the field values at points A, B, E, and F in Fig. £|. The first two of these 
field values are already known at the nth step, while the last two can be obtained by cubic 
radial interpolations along the previously evolved characteristics u = n ra _i and u = n n _ 2 , 
respectively. At the initial step, point F lies on the characteristic u = n_i = —At — R m , 
which is not evolved by the algorithm; this field value is supplied along with the initial data. 
Once the Cauchy and characteristic boundary values are computed, a new iteration may be 
performed starting from Step 1 above. 

Since all the interpolations employed in the matching step have fourth order error, the 
matching algorithm has the same second order global accuracy exhibited by the separate 
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Cauchy and characteristic algorithms, as confirmed by numerical tests. 

In summary, CCM has been implemented and tested for nonlinear waves in 3-dimensional 
space. No special assumption is made about the waves crossing the computational interface 
r = R m and nonlinear effects in the exterior characteristic domain are automatically taken 
into account. In numerical experiments CCM converged to the exact solution (with the 
matching boundary fixed at an arbitrary position) in highly nonlinear problems. For com- 
parison, nonlocal ABC's yielded convergent results only in linear problems. In terms of both 
computational cost and accuracy, CCM is a very effective way to solve the 3-dimensional 
wave equations, with or without a nonlinear term. It is possible to achieve convergence with 
an ABC by refining the grid and simultaneously enlarging the radius of the outer bound- 
ary. However, this is very expensive computationally, especially for small target error in 
the determination of the radiated signal in a 3-dimensional problems [[31]. Because CCM is 
convergent under grid-refinement alone, for small target error its performance is significantly 
better than any available alternative. In strongly nonlinear problems. CCM appears to be 
the only available method which is able to produce numerical solutions which converge to 
the exact solution with a computational interface located at an arbitrary fixed position. 



V. APPLICATION TO BLACK HOLES 

By definition, a black hole traces out a world tube in space-time which is the boundary 
of what is visible to an outside observer. As a result, this worldtube is necessarily a char- 
acteristic hypersurface traced out in space-time by light rays (the characteristics of general 
relativity). 

A mechanical analogue of a black hole can in principle be constructed by letting a 
reservoir of water empty through a funnel into a lower reservoir. By arranging the flow 
velocity in the funnel to exceed the velocity of sound in the upper reservoir, an acoustic 



version of a black hole results ||32|| . Sound waves can travel from the upper reservoir to the 
lower but not in the reverse direction. 

The simplest example of a black hole arises in the spherically symmetric collapse of a 
star whose energy has been depleted to the extent that internal pressure cannot withstand 
gravity. A point sized black hole first forms at the center of the collapsing star, which then 
expands into a sphere growing with the velocity of light, thus tracing out a light cone in 
space-time. Normally a light cone keeps expanding forever. What makes a black hole light 
cone unique is that gravity halts this expansion and eventually the black hole just hovers in 
equilibrium at a fixed size. Each point on the spherical black hole still moves along a light 
ray, but the sphere of light rays is in a delicate gravitational balance between growing and 
shrinking. 

This spherically symmetric black hole was discovered analytically by Schwarzschild as a 
simple solution to Einstein's equations. Unfortunately, Schwarzschild's black hole does not 
emit gravitational waves (just as a spherically symmetric charge distribution does not emit 
electromagnetic waves). 

The inspiral and merger of a binary system of black holes is a powerful source of grav- 
itational waves, which lie in the frequency band detectable by the new gravity wave obser- 
vatories. The lack of symmetry necessitates a computational treatment to determine the 
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waveform of the radiated signal. The binary black holes trace out a characteristic hypersur- 
face in space-time so that their simulation by characteristic evolution is a natural approach. 
Using a characteristic evolution algorithm similar to that which I have described here, ex- 
cept now applied to Einstein's equation, we have found fascinating results. The binary black 
holes, rather than initially forming as a point caustic in the Schwarzschild case, form as a 
cross-over surface where light rays collide. This cross-over surface is itself bounded by a ring 
of caustics which in a sense mark the merger of the two individual black holes into a single 



black hole. 33.34 



The individual black holes form as spheres but, as illustrated in Fig. ^j, as the holes 
approach, their mutual gravitational tidal distortion produces sharp pincers just prior to 
merger. At merger, the pincers join to form a single temporarily toroidal black hole, as 
illustrated in Fig. |. The inner hole of the torus subsequently closes up to produce first a 
single peanut shaped black hole and finally a spherical black hole. Details of this merger 
can be viewed at the web site [http: / / artemis.phyast.pitt.edu/ animations!- We are now using 



characteristic evolution to compute the gravitational wave signal emitted in such black hole 
collisions in the anticipation that they will be detected by the new gravity wave observatories. 
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FIGURES 



u 




FIG. 1. Space-time picture of a parallelogram E formed by intersecting characteristics, where 
the characteristic have 45° slope. The coordinates are retarded time u and radial distance r. In 
standard (t, x) coordinates, the t-direction is vertical and the spatial x-direction horizontal, so that 
the pair P and S are at the same time t, with Q in their future and S in their past. In (u, r) 
coordinates, the pair P and Q are at the same retarded time, as well as the pair R and S. 
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FIG. 2. Line segments drawn at forty-five degrees represent radial characteristics. Their inter- 
section defines the fundamental null parallelogram PQRS shown superimposed upon the computa- 
tional cell, which consists of the points marked by circles and their nearest neighbors in the angular 
directions (not shown). Here u is the retarded time coordinate and x is the compactified radial 
coordinate. The marching algorithm determines <£q in terms of the previously determined values 
3>P, <&r and $5. The process is then iterated to determine <&t on a march to radiative infinity. 
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FIG. 3. Initial Cauchy data are evolved from to to time t\ throughout the region D]_- . Char- 
acteristic data induced on C x - , combined with the initial characteristic data on C + are used to 
evolve the region D^. This produces Cauchy data at time t\ in the region r < R m . Similarly, 
Cauchy evolution is used in the region D 2 - , bounded on the right by C 2 - . The characteristic data 
induced on C 2 - , together with those on C\+ , are sufficient to evolve through the region D 2 + ■ The 
process can be iterated to carry out the entire future evolution of the system. 
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FIG. 4. Cauchy grid points are indicated by squares and characteristic grid points by circles. 
The triangles indicate points E and F where the time level t n intersects the retarded time levels 
-u n _i and u n -2- Initial data are given at the shaded points. Evolution proceeds iteratively by 
determining field values at the unshaded points. The matching scheme provides boundary values 
at C (r = R m ) and D (r = Rb + Ar) for the characteristic and Cauchy grids, respectively. 
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FIG. 5. Prelude to a generic collision of black holes. As they approach, the initially spherical 
black holes are tidally distorted. 




FIG. 6. Temporarily toroidal black hole produced by the merger of two black holes. The hole 
in the torus soon closes up to form a peanut shaped black hole, which then expands into a spherical 
shape. 
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